1 Contexto

O cancelamento de assinatura é uma pedra no sapato de empresas de telecomunicações e streaming, por exemplo, pois traz instabilidade financeira e prejudica o planejamento estratégico da companhia. Por isso, fazer um acompanhamento da taxa de cancelamento ou Churn Rate é fundamental para manter o negócio vivo e de pé. Mas além de monitorar o indicador, é possível prever quando um cliente deve cancelar sua assinatura? Sim, é possível e é isso que vamos mostrar logo mais. Seja bem vindo e te desejo uma boa leitura!

1.1 Origem dos dados

Os dados que vamos utilizar para aplicar aprendizado de máquina estão aqui.

1.2 Objetivo

Esta análise tem por finalidade prever o cancelamento de assinatura por parte do cliente, aplicando um modelo de aprendizado de máquina.

2 Pacotes

library(DescTools)
library(tidyverse)
library(tidymodels)
library(reactable)
library(janitor)
library(highcharter)
library(sparkline)
library(corrplot)
library(patchwork)
library(rmarkdown)
library(knitr)
library(kableExtra)
library(upstartr)
library(tictoc)
library(vip)
theme_set(theme_light(18))

3 Carregando os dados

dados = tibble(read_csv("C:/Users/icaro/Documents/Meu_Portfolio/Churn_Telecom/WA_Fn-UseC_-Telco-Customer-Churn.csv"))

O que temos dentro desse conjunto?

paged_table(dados)

Então temos um dataset com 7043 registros/observações e 21 atributos ou melhor, variáveis.

4 Exploratória (EDA)

Vamos explorar nosso dataset, começando por conhecer nossas variáveis.

# contagem
vars = 
  tibble(
    tipos = c("quantitativa", "qualitativa"),
    colunas = 
      c(
        dados %>% select(where(is.numeric)) %>% ncol(),
        dados %>% select(where(is.character)) %>% ncol()
    ))

# porcentagem
vars = 
  vars %>% 
  mutate(
    porcent = round(colunas / sum(colunas), digits = 2)
  )

# rótulos
rotulo = paste(str_to_title(vars$tipos), " (", vars$porcent * 100, "%", ")", sep = "")

# gráfico
pie(
  x = vars$colunas, 
  labels = rotulo, 
  col = c("skyblue2", "skyblue4"),
  border = "white"
)
title(
  main = "Tipo de Variáveis", 
  cex.main = 2, 
  font.main = 4, 
  col.main = "skyblue4"
  )

Temos que em sua maior parte, nosso conjunto de dados é formado por variável qualitativa nominal.

Antes de iniciarmos a análise exploratória, vamos realizar alguns ajustes no dataset.

Excluir o campo de ID e padronizar o nome das variáveis.
Transformando o campo Senior Citizen em string.

# padronizando os nomes da vars, excluir a coluna com o ID
dados = 
  dados %>% select(-customerID) %>% clean_names()

# formatando o campo senior citizen
dados$senior_citizen = ifelse(dados$senior_citizen == 0, 'No', 'Yes')

4.1 Sumário

Agora sim, vamos realizar a construção de um sumário separando entre variáveis numérica e categórica.

# Sumário - Variáveis Categóricas
sumario_character = 
  dados %>%
  select(where(is.character)) %>%
  as.list() %>%
  enframe(name = "variavel", value = "valores") %>%
  mutate(
    n_missing = map_dbl(valores, ~sum(is.na(.x))),
    complete_rate = 1 - n_missing/map_dbl(valores, ~length(.x)),
    min = map_dbl(valores, ~min(str_length(.x))),
    max = map_dbl(valores, ~max(str_length(.x))),
    # empty = map_dbl(valores, ~sum(str_detect("^$", .x))),
    n_unique = map_dbl(valores, ~n_distinct(.x)),
    whitespace = map_dbl(valores, ~sum(str_detect("^[:blank:]+$", .x)))
  ) %>%
  arrange(variavel != "churn") 

sumario_character %>%
  reactable(
    wrap = FALSE, 
    resizable = TRUE,
    fullWidth = TRUE,
    defaultColDef = colDef(width = 90),
    columns = list(
      valores = colDef(show = FALSE),
      variavel = colDef("Variável", minWidth = 230, width = 230)
    ),
    details = function(index) {
      variavel_chr <- sumario_character[index, "variavel", drop = TRUE]
      dados %>%
        tabyl(!!sym(variavel_chr)) %>% arrange(desc(n)) %>%
        reactable(columns = list(percent = colDef("%", format = colFormat(percent = TRUE, digits = 1))), width = 500)
    }
  )

Agora para as variáveis numéricas.

# Sumário - Variáveis Numéricas
sumario_numeric <- dados %>%
  select(where(is.numeric)) %>%
  as.list() %>%
  enframe(name = "variavel", value = "valores") %>%
  mutate(
    n_missing = map_dbl(valores, ~sum(is.na(.x))),
    complete_rate = 1 - n_missing/map_dbl(valores, ~length(.x)),
    mean = map_dbl(valores, ~mean(.x)),
    sd = map_dbl(valores, ~sd(.x)),
    min = map_dbl(valores, ~min(.x)),
    median = map_dbl(valores, ~median(.x)),
    max = map_dbl(valores, ~max(.x))
  ) %>%
  mutate(across(where(is.numeric), round, digits = 1))

sumario_numeric %>%
  reactable(
    wrap = FALSE, 
    resizable = TRUE,
    fullWidth = TRUE,
    defaultColDef = colDef(width = 90),
    columns = list(
      valores = colDef(cell = function(values) {sparkline(table(cut(values, min(n_distinct(values), 15))), type = "bar")}),
      variavel = colDef("Variável", minWidth = 230, width = 230)
    ),
    details = function(index) {
      hchart(sumario_numeric[index, "valores"][[1]][[1]])
    }
  )

4.2 Relação entre as explicativas

4.2.1 Correlação variáveis numéricas

# Correlações entre as numéricas
dados %>% 
  select(where(is.numeric)) %>% 
  filter(!is.na(dados$total_charges)) %>%  
  cor() %>% 
  corrplot(method = "color", 
           order = "hclust", 
           type = "upper", 
           addCoef.col = "black",
           tl.col = "black",
           cl.align.text = "l")

4.2.2 Sankey explicativas categóricas

# Sankey das explicativas categóricas
dados %>%
  select(where(is.character), -churn) %>%
  data_to_sankey() %>%
  hchart("sankey")

4.2.3 Relação com a variável resposta

# base transformada
churn_numeric_long = 
  dados %>%
  select(where(is.numeric), churn) %>%
  pivot_longer(-churn, names_to = "variavel", values_to = "valor")

4.2.3.1 Variáveis numéricas

# plot base
plot_base = 
  churn_numeric_long %>%
  ggplot(aes(x = valor)) +
  facet_wrap(~variavel, scale = "free", ncol = 1) +
  theme(axis.text = element_blank(),
        axis.title = element_blank())

# Histogramas
p1 <- plot_base + geom_histogram(aes(fill = churn), position = "identity")

# Densidades
p2  <- plot_base + geom_density(aes(fill = churn), alpha = 0.4, show.legend = FALSE)

# Boxplots
p3 <- plot_base + geom_boxplot(aes(fill = churn), show.legend = FALSE)

# KS
p4 <- plot_base + stat_ecdf(aes(colour = churn), show.legend = FALSE)

# plot all
(p4 + p3 + p2 + p1) + plot_layout(nrow = 1)

4.2.3.2 Variáveis categóricas

churn_character_long = 
  dados %>%
  select(where(is.character), churn) %>%
  pivot_longer(-churn, names_to = "variavel", values_to = "valor")

# bar plot base
bar_plot = 
  churn_character_long %>%
  count(variavel, valor, churn) %>%
  group_by(variavel, valor) %>%
  mutate(
    p_churn = sum(n[churn == "sim"], na.rm = TRUE)/sum(n, na.rm = TRUE)
  ) %>%
  group_by(variavel) %>%
  mutate(
    valor = reorder_within(valor, p_churn, variavel)
  ) %>%
  ggplot(aes(y = valor, x = n, fill = churn)) +
  facet_wrap(~variavel, scales = "free", ncol = 1) +
  scale_y_reordered() +
  theme(axis.text.x = element_blank(),
        axis.title = element_blank())

# Position stack
p1 = bar_plot + geom_col(show.legend = FALSE, position = "stack")

# Position fill
p2 = bar_plot + geom_col(show.legend = TRUE, position = "fill")

p1 + p2

5 Transformação

Vimos que há variáveis categóricas que precisam de alteração para fator, além agrupar os casos de clientes que não possuem serviço de internet ao grupo dos ‘NO’, pois representam a “mesma” situação.

# tranformando em fator as variaveis categoricas
dados_transf = 
  dados %>% 
  mutate(
    across(
      c(gender:dependents, 
        phone_service, multiple_lines,
        online_security:streaming_movies,
        paperless_billing, churn), 
      factor),
    multiple_lines = 
      ifelse(multiple_lines == 'No phone service', 'No', 
             ifelse(multiple_lines == 'No', 'No', 'Yes')),
    online_security = 
      ifelse(online_security == 'No internet service', 'No', 
             ifelse(online_security == 'No', 'No', 'Yes')),
    online_backup = 
      ifelse(online_backup == 'No internet service', 'No', 
             ifelse(online_backup == 'No', 'No', 'Yes')),
    device_protection = 
      ifelse(device_protection == 'No internet service', 'No', 
             ifelse(device_protection == 'No', 'No', 'Yes')),
    tech_support = 
      ifelse(tech_support == 'No internet service', 'No', 
             ifelse(tech_support == 'No', 'No', 'Yes')),
    streaming_tv = 
      ifelse(streaming_tv == 'No internet service', 'No', 
             ifelse(streaming_tv == 'No', 'No', 'Yes')),
    streaming_movies = 
      ifelse(streaming_movies == 'No internet service', 'No', 
             ifelse(streaming_movies == 'No', 'No', 'Yes'))
    )

6 Modelagem

Vamos iniciar o processo de modelagem dos dados com a separação da base em teste e treino.

# marcando a base em teste e treino
set.seed(1232)
churn_split = 
  initial_split(data = dados_transf, strata = churn, prop = 0.75)

# separando a base em treino
training_df = training(churn_split)
# separando a base em teste
testing_df = testing(churn_split)

6.1 Validação cruzada

# avaliação dos modelos com cross-validation
folds = 
  vfold_cv(training_df, v = 10, strata = churn)

6.2 Receita

Na receita informamos as etapas que precisamos realizar para preparar a base de treino do modelo.

No nosso caso, a receita inclui:

Correção do balanceamento de classes.
Exclusão de variáveis com alto coeficiente de correlação.
Normalização das variáveis numéricas.
Substituição de valores faltantes pela mediana em campos númericos.
Transformação das variáveis qualitativas nominais em dummy.

# preparando os dados 
churn_rec = 
  recipe(churn ~ ., data = training_df) %>% 
  themis::step_downsample(churn) %>% 
  step_corr(all_numeric_predictors()) %>% # correlação entre explicativas
  step_normalize(all_numeric()) %>% # normalizar vars
  step_rm(total_charges) %>% # Tratamento de missings
  step_impute_median(all_numeric()) %>% # Substituir os missings com a mediana 
  step_dummy(all_nominal(), -churn) # Transformando as variaveis nominais em dummy
  
# base processada  
baked = bake(prep(churn_rec), new_data = NULL)
paged_table(baked)

6.3 Modelo Random Forest

# modelo de randow forest (árvore aleatória)
rf_mod = 
  rand_forest(
    trees = 300,
    mtry = tune(),
    min_n = tune()
  ) %>%
  set_mode('classification') %>% 
  set_engine('ranger', importance = 'impurity')

6.4 Workflow

rf_workflow = 
  workflow() %>% 
  add_recipe(churn_rec) %>% 
  add_model(rf_mod)

6.5 Tunagem

set.seed(1)
tic("modelo rf")
tunagem = 
  tune_grid(
    rf_workflow,
    resamples = folds,
    grid = 10,
    metrics = metric_set(roc_auc, precision, accuracy, f_meas),
    control = control_grid(verbose = TRUE, allow_par = FALSE)
    )
toc()

## modelo rf: 79.81 sec elapsed

6.6 Avaliações

# tabela
show_best(tunagem, "roc_auc") %>% 
  kable(format = 'html', align = 'l', digits = 3) %>% 
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

mtry	min_n	.metric	.estimator	mean	n	std_err	.config
8	31	roc_auc	binary	0.848	10	0.005	Preprocessor1_Model02
13	36	roc_auc	binary	0.848	10	0.005	Preprocessor1_Model01
22	38	roc_auc	binary	0.845	10	0.005	Preprocessor1_Model10
4	8	roc_auc	binary	0.844	10	0.006	Preprocessor1_Model09
17	26	roc_auc	binary	0.844	10	0.005	Preprocessor1_Model04

6.7 Atualização do workflow

Atualizando o workflow com os hiperparamêtros encontrados.

# atualizando wf
wf = 
  rf_workflow %>% 
  finalize_workflow(select_best(tunagem, "roc_auc"))

6.8 Ajuste

# last fit
ajuste_final = 
  last_fit(
    wf, 
    churn_split, 
    metrics = 
      metric_set(accuracy, roc_auc, f_meas, specificity, precision, recall)
    )

# métricas de desempenho
collect_metrics(ajuste_final) %>% 
  kable(format = 'html', digits = 3, align = 'l') %>% 
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

.metric	.estimator	.estimate	.config
accuracy	binary	0.740	Preprocessor1_Model1
f_meas	binary	0.805	Preprocessor1_Model1
specificity	binary	0.769	Preprocessor1_Model1
precision	binary	0.897	Preprocessor1_Model1
recall	binary	0.730	Preprocessor1_Model1
roc_auc	binary	0.825	Preprocessor1_Model1

# predicoes
pred_base_teste = collect_predictions(ajuste_final)

# curva roc
roc = 
  pred_base_teste %>% 
  roc_curve(churn, .pred_No) %>%
  autoplot()

# curva de lift
lift = 
  pred_base_teste %>% 
  lift_curve(churn, .pred_No) %>%
  autoplot()

# KS
ks = 
  pred_base_teste %>% 
  ggplot(aes(x = .pred_Yes, colour = churn)) +
  stat_ecdf(show.legend = FALSE)
  
# distribuicao
dist = 
  pred_base_teste %>% 
  ggplot(aes(x = .pred_Yes, fill = churn)) +
  geom_density() +
  theme(axis.title = element_blank())
  
(roc + lift)/(ks + dist)

A conclusão que tiramos daqui é que o modelo apresentou razoável acurácia, separou bem os grupos, mas há espaço para melhorar, pois existe um ponto central no gráfico de densidade que mostra que o modelo não conseguiu acertar.

6.9 Modelo final

Ajuste final com a base inteira.

modelo_final = fit(wf, dados_transf)

6.10 Importância das variáveis

Avaliação das contribuição das variáveis para o modelo.

vip(modelo_final$fit$fit)

6.11 Armazenamento

Podemos salvar o modelo para realizar as predições futuramente, sem a necessidade de rodar o script sempre que precisar.

saveRDS(modelo_final, file = "modelo_final.rds")

6.12 Matrix de confusão

Vamos confrontar as previsões do modelo com a realidade.

predizer = 
  tibble(
    cbind(
      dados_transf$churn, 
      predict(modelo_final, new_data = dados_transf, type = 'prob')
      )
  )

predizer = 
  predizer %>% 
  rename(churn_observado = `dados_transf$churn`) %>% 
  mutate(churn_esperado = ifelse(.pred_No > .pred_Yes, 'No', 'Yes'),
         .pred_No = round(.pred_No, digits = 3),
         .pred_Yes = round(.pred_Yes, digits = 3))

# matriz
predizer %>% 
  tabyl(churn_observado, churn_esperado) %>% 
  adorn_percentages(denominator = "all") %>%  
  adorn_pct_formatting() %>% 
  adorn_title(row_name = "Observado", col_name = "Esperado") %>% 
  kable(format = "html", align = 'c', caption = "Matriz de confusão") %>% 
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Matriz de confusão
	Esperado
Observado	No	Yes
No	56.9%	16.6%
Yes	3.5%	23.0%

6.13 Predição

Vamos testar a predição do modelo, inserindo os dados de um ‘novo cliente’.

novo_cliente = 
  data.frame(
    gender = 'Female',
    senior_citizen = 'Yes',
    partner = 'Yes',
    dependents = 'Yes',
    tenure = 25,
    phone_service = 'No',
    multiple_lines = 'No',
    internet_service = 'DSL',
    online_security = 'No',
    online_backup = 'No',
    device_protection = 'No',
    tech_support = 'No',
    streaming_tv = 'Yes',
    streaming_movies = 'Yes',
    contract = 'Month-to-month',
    paperless_billing = 'Yes',
    payment_method = 'Electronic check',
    monthly_charges = 55.55,
    total_charges = 290.50
    )


predict(modelo_final, new_data = novo_cliente, type = 'prob') %>% 
  kable(
    format = "html", 
    digits = 3, 
    align = "c", 
    caption = "Novo cliente",
    col.names = c("Não", "Sim")
    ) %>% 
  kable_styling(
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover", "condensed", "responsive")
    )

Novo cliente
Não	Sim
0.509	0.491

Neste exemplo, com a probabilidade maior para não, então assumimos que o cliente não deve cancelar sua assinatura.

7 Final

Recaptulando o fizemos aqui até.

Análise descritiva da base, explorando as relações entre as variáveis.
Transformação necessária para aplicação do modelo.
Modelagem, com a especificação da f(x) - receita - workflow até o armazenamento do modelo e avaliação.

Ao final deste estudo, podemos concluir que o modelo conseguiu realizar um trabalho razoável, com precisão de 0,89. No entanto, a matriz de confusão nos mostrou que o modelo errou bastante, onde previu que o cliente cancelaria a assinatura quando na verdade não cancelou.

Churn Telecom

Predizendo o cancelamento de assinatura

Ícaro Pinheiro - Estatístico e Analista de Dados

2022-11-12